/*=====================================================================*
 *                   Copyright (C) 2011 Paul Mineiro                   *
 * All rights reserved.                                                *
 *                                                                     *
 * Redistribution and use in source and binary forms, with             *
 * or without modification, are permitted provided that the            *
 * following conditions are met:                                       *
 *                                                                     *
 *     * Redistributions of source code must retain the                *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer.                                       *
 *                                                                     *
 *     * Redistributions in binary form must reproduce the             *
 *     above copyright notice, this list of conditions and             *
 *     the following disclaimer in the documentation and/or            *
 *     other materials provided with the distribution.                 *
 *                                                                     *
 *     * Neither the name of Paul Mineiro nor the names                *
 *     of other contributors may be used to endorse or promote         *
 *     products derived from this software without specific            *
 *     prior written permission.                                       *
 *                                                                     *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND              *
 * CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES,         *
 * INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES               *
 * OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE             *
 * ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER               *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,                 *
 * INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES            *
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE           *
 * GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR                *
 * BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF          *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT           *
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY              *
 * OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE             *
 * POSSIBILITY OF SUCH DAMAGE.                                         *
 *                                                                     *
 * Contact: Paul Mineiro <paul@mineiro.com>                            *
 *=====================================================================*/

#include "wexpxmx.hh"

namespace
{
  float
  newton (float x,
          float z)
    {
      // z should be a solution to

      // x + z - Exp[-z] == 0

      float expz = expf (z);

      return (1 - expz * x + z) / (1 + expz);
    }

  float
  poly (float        x,
        float*       coeffs,
        unsigned int n)
    {
      float z = coeffs[0];

      for (unsigned int i = 1; i <= n; ++i)
        {
          z *= x;
          z += coeffs[i];
        }

      return z;
    }

  float 
  wexpxmx_m8_m6 (float x) 
    {
      float coeffs[] = { -7.326416836769547e-7,
                         -0.000021982630104104322,
                         -0.000247896530953556,
                         -1.001245857586513,
                         -0.00235587292524 };

      return expf (x) + poly (x, coeffs, 4);
    }

  float 
  wexpxmx_m6_m2 (float x) 
    {
      float xplusfour = x + 4.0;
      float numerator[] = { 9.8587764254740447e-7,
                            0.000015156527458371576,
                            0.00013331812995983005,
                            -0.025484803751509787,
                            0.31953133142385250,
                            -1.6510673002491159,
                            3.2938619987374864 };
      float denominator[] = { 0.026322737102276157,
                              -0.21049708859450385,
			      0.81977872872246388 };
                         
      float z = 
        poly (xplusfour, numerator, 6) / poly (xplusfour, denominator, 2);

      if (x > -2.75)
        {
          z = newton (x, z);
        }

      return z;
    }

  float
  wexpxmx_m2_m0 (float x)
    {
      float xplusone = x + 1.0;
      float numerator[] = { -0.000020329490994144038,
                            0.000079465027615349629,
                            0.0027610311594243719,
                            -0.045474624127473696,
                            0.15440546182117503,
                            -0.84315650114928865,
                            1.3870842152670159 };

      float denominator[] = { 0.066854477487842454,
                              0.0042921144767478142,
                              1.0849610363626968 };

      return poly (xplusone, numerator, 6) / poly (xplusone, denominator, 2);
    }

  float
  wexpxmx_0_2 (float x)
    {
      float xminusone = x - 1.0;
      float numerator[] = { 4.6422073954430867e-6,
                            -0.000099888341204264827,
                            0.0012187632506108073,
                            -0.012590998049710100,
                            -0.058701377425761985,
                            -0.50753166974511571,
                            0 };
      float denominator[] = { 0.045144135953763134,
                              0.24428567225939440,
                              1.0150633392629635 };

      return poly (xminusone, numerator, 6) / poly (xminusone, denominator, 2);
    }

  float
  wexpxmx_2_4 (float x)
    {
      float xminusthree = x - 3.0;

      float numerator[] = { 6.2024884510996672e-7,
                            -0.000012192552084561388,
                            0.00016403486040872973,
                            -0.0023513476270810808,
                            -0.054959349893915783,
                            -0.47679270789600670,
                            -0.75561412175512671 };
      float denominator[] = { 0.020518884102674959,
                              0.22651058466352914,
                              0.95398600089869317 };

      return 
        poly (xminusthree, numerator, 6) / poly (xminusthree, denominator, 2);
    }

  float
  wexpxmx_4_8 (float x)
    {
      float xminussix = x - 6.0;

      float numerator[] = { 1.3776323249096728e-8,
                            -4.6623963046723786e-7,
                            0.000012424071435159243,
                            -0.00043990787481529912,
                            -0.040247069742738717,
                            -0.47953110077434266,
                            -1.4560399693642404 };

      float denominator[] = { 0.011076647304105349,
                              0.20176872873568395,
                              0.96853939300829810 };

      return poly (xminussix, numerator, 6) / poly (xminussix, denominator, 2);
    }

  float
  wexpxmx_gt8 (float x)
    {
      float logx = logf (x);

      float z = -logx + logx / x + (logx - 2) * logx / (2 * x * x);

      if (x < 20)
        {
          z = newton (x, z);
        }

      return newton (x, z);
    }
};

namespace flassol
{
  float 
  wexpxmx (float x)
    {
      return (x < -16) ? -x :
             (x < -8) ? -x + expf (x) :
             (x < -6) ? wexpxmx_m8_m6 (x) :
             (x < -2) ? wexpxmx_m6_m2 (x) :
             (x < 0) ? wexpxmx_m2_m0 (x) :
             (x < 2) ? wexpxmx_0_2 (x) :
             (x < 4) ? wexpxmx_2_4 (x) :
             (x < 8) ? wexpxmx_4_8 (x) :
                       wexpxmx_gt8 (x);
    }
}


